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ABSTRACT 


The Helmholtz Equation 
(-A-K^n^)u = 0 

with a variable index of refraction, n, and a suitable radiation condition at 
infinity serves as a model for a wide variety of wave propagation problems. A 
numerical algorithm has been developed and a computer code Implemented that 
can effectively solve this equation in the intermediate frequency range. The 
equation is discretized using the finite element method, thus allowing for the 
modeling of complicated geometries (including interfaces) and complicated 
boundary conditions. A global radiation boundary condition is Imposed at the 
far field boundary that is exact for an arbitrary number of propagating modes. 

The resulting large, non-self adjoint system of linear equations with 
indefinite symmetric part is solved using the preconditioned conjugate 
gradient method applied to the normal equations. A new preconditioner is 
developed based on the multigrid method. This preconditioner is vectorizable 
and is extremely effective over a wide range of frequencies provided the 
number of grid levels is reduced for large frequencies. A heuristic argument 
is given that Indicates the superior convergence properties of this 
preconditioner. The relevant limit to analyze convergence is for K 
increasing and a fixed prescribed accuracy level. The efficiency and 
robustness of the numerical algorithm is confirmed for large acoustic models, 
including Interfaces with strong velocity contrasts. 
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Introduction 


In this paper, we describe a numerical method for approximately solving 
the Helmholtz equation 

Au + K^n^(x,y)u = 0. (1»1) 

Equation (1.1) with a suitable radiation condition at infinity describes both 

the propagation and scattering of time harmonic waves in general geometries. 

We will restrict the application of (1.1) to problems that occur in underwater 

acoustics. Therefore u will represent the acoustic pressure, n(x,y) the 

2tt f 

index of refraction, and K is the wave number (- , where f is the 

frequency and Cq is a reference sound speed). The region of Interest will 
be a duct or waveguide containing inhomogeneities and interfaces. The 
numerical method is based on combining a finite element discretization with a 
recently developed iterative method for solving the resulting system of linear 
equations [1]. We observe that this numerical method is also applicable to 
three-dimensional problems as well as vector formulations of (1.1) such as 
those describing elastic wave propagation [2]. 

Various computational techniques have been applied to simplified 
propagation models, including parabolic equation and normal mode type methods, 
asymptotic methods, and others. For a survey of various models and 
computational methods, see [3]. Although each of these techniques is 

effective under suitable assumptions, there are many important problems for 
which it is necessary to treat the complete wave propagation model in the low 
to intermediate frequency range. Such models can include full angle 
propagation and backscatterlng. This can occur, for example, when the ocean 



- 2 - 


bottom must be taken Into account or a layer of ice is present on the ocean 
surface. 

There are several difficulties associated with the numerical solution of 
the full acoustic propagation model (1.1). The solution must often be 
resolved over many wavelengths, leading to very large systems of linear 
equations. Recent results have shown that as the frequency Increases (or 
equivalently as the size of the model Increases), the number of 

points/wavelength must be Increased [4]. Thus the number of equations 
increases faster than quadratlcally in K. A finite element method has 
recently been developed that dramatically reduces the number of equations in 
regions where little backscatterlng is present [5]. The resulting matrices 
will be indefinite and also nonself-adjoint due to the radiation boundary 
condition. Furthermore, effective radiation boundary conditions to be imposed 
at a finite boundary must be developed. Various alternative approaches for 
dealing with these difficulties are discussed in [5]. 

We have developed a numerical algorithm and Implemented a computer code 
that can effectively solve (1.1). The equation is discretized using the 
finite element method, thus allowing considerable flexibility in modeling 
complex geometries. The resulting linear system of equations is solved 
iteratively using the preconditioned conjugate gradient method applied to the 
normal equations [1]. This method requires relatively little storage' (l.e. , 
storage does not have to be allocated for the bandwidth) and is well suited 
and efficient for large problems. A global radiation boundary condition is 
Imposed at the far field boundary. This boundary condition is based on a 
modal expansion of the far field solution that is exact for an arbitrary 
number of propagating modes (see [6] and [7]). 


- 3 - 


The effectiveness of the iterative method depends on the choice of the 
preconditioner. We consider preconditioners based on a splitting of the 
discrete Laplacian. In previous papers, [1] and [8], we investigated 
preconditioners based on SSOR and ADI. In this paper we describe a 
preconditioner based on a version of the multigrid method (introduced in [16]) 
which we have found to be extremely effective over a wide range of 
frequencies, provided the number of grid levels is decreased as the frequency 
increases. We shall see that this preconditioning is considerably more 
efficient than the SSOR and ADI preconditioners employed in [1] and [8]. This 
preconditioning has the additional advantage of being vectorizable since a 
relaxation scheme based on a red-black ordering [9] is used. 

We close this section by outlining the remainder of the paper. In Section 
2 we describe a wave propagation model Including an interface with a 
discontinuity in both the sound speed and density. We also describe the 
global radiation boundary condition, a finite element method for solving this 
boundary value problem, and an efficient implementation of the global boundary 
condition condition with the finite element method. In Section 3 we describe 
the iterative solution method and multigrid preconditioner. We also describe 
a new way of analyzing convergence of the iterative method, based on an 
accuracy condition developed in [4] relating the frequency and mesh size. 
Numerical results will be presented in Section 4 demonstrating the efficiency 
of the numerical method. We summarize our conclusions in Section 5. 
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Contlnuous and Discrete Model 

We shall describe our numerical method with respect to the following model 
problem. Consider (1.1) in a bounded portion of a two-dimensional semi- 
infinite rectangular waveguide. Hence our computational domain, D, is given 
by D = {(x,y); 0 < x < x^, 0 < y < ir} . We assume that there is an Interface 
r dividing D into two subregions, and D 2 . Furthermore, suppose that the 
density p is piecewise constant with p = and and p = P 2 in D 2 . 

Our propagation model is now given by the following boundary value problem: 

(-A - K^n^(x,y))u(x,y) = 0 in D (2.1a) 

3u(x,0) ^ 0 (2.1b) 

dy 

u(x,it) = 0 (2.1c) 


9u(0,y) 

3x 


= g(y) 


(2. Id) 


. T(u) 


(2.1e) 


Uj(x,y) = U 2 (x,y), for (x,y) on T 


(2. If) 


and 

, 3u,(x,y) , 3u2(x,y) 

P 2 for (x,y) on T, (2.1g) 


The boundary operator T in (2.1e) is chosen to model the outflow of energy 
and will be discussed in detail below. In (2. If) and (2.1g), u^ denotes the 
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restriction of the pressure u to D^, i=l,2, and 3uj^/3n denotes the normal 
derivative on T pointing into ^2* 


Remark 2.1 . The Dirichlet boundary condition (2.1c) is a pressure release 
condition valid on the ocean surface. Condition (2.1b) models a rigid bottom, 
although a more general Impedance condition could be Implemented without 
difficulty. The forcing term in (2. Id) could readily be imposed as a 
Dirichlet condition Instead of a Neumann condition. Furthermore, we could 
just as easily consider other coordinates systems (e.g., cylindrical 
coordinates) instead of Cartesian coordinates. Finally, note that the method 
can be applied to problems with nonrectangular boundaries (see [6]). Based on 
our previous experience with complicated geometries [8], we do not anticipate 
a serious degradation in the numerical results. 

We next define the radiation boundary operator T appearing in (2.1e). 
Consider the semi-infinite rectangular waveguide 

W = {(x,y): 0<x<«», 0<y<w}. Assuming that n(x,y) = I for x > x , 
it is easily seen that the outgoing solution of (2.1) can be expressed as 


u(x,y) 



> x^, where = j +V 2 and 

for K > Oj 
for < Oj 


Thus the far field solution is composed of a finite number of progagatlng 
modes (Ij imaginary) and an infinite number of evanescent modes (i^ real) 
that decay exponentially as x Set f .(y) = cos o.y and note that 

1 j 


T 
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a. = e ^ “ < u,f > = e J “ / u(x^,y)f (y)dy, j=l,2,»»«, 

J J P 0 

00 

where 

= {(x,y): X = , 0 < y < it}. 

We may now define our global boudnary operator, as in [6] and [7], by 

T(u) = I <«>Vr 
m=l « 

where M is chosen large enough to account for all propagating inodes and 
those evanescent modes which are not small enough to be neglected at x = x^. 

Remark 2 » 2 : When only very few propagating modes are important, a simpler 

local boundary operator can be efficiently implemented [17, 18] • For example, 

£ X 

if only the mth mode, e f^(y), is propagating for x > x^, the global 
operator defined by (2.2) could be replaced by the impedance boundary 

operator given by 

B(u) = (|^ - iA^)u. (2.2') 

However, as the number of significant modes increases, the order of the local 
boundary operator Increases. 

We next discretize the continuous model, (2.1) using the finite element 
method. We first observe, using integration by parts, that the solution u 

of (2.1) satisfies the following variational problem: 

IT 

a(u,v) = -p_<g,v>„ = -p- / g(0,y)v(0,y) dy for all v in (2.3) 

0 0 
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(assuming C D^, where s {(x,y); x = 0, 0 < y < ir} and is the 

space of continuous complex-valued functions defined by 


H 


E _ 


{v: 


2 

I 


/ (M' 


2 

+ |Vv| )dxdy < • and v(x,ir) - o}. 


The bilinear form a(v,w) is defined by 


a(v,w) = p„ / (Vv»Vw - K n w)dxdy 


D. 


+ Pj(/ (Vv*Vw - K^n^vw)dxdy - ^ T(v)wdy) 


(2.4) 


for all v,w in H 


E 


To discretize problem (2.3) we introduce a family of finite dimensional 
subspaces such that becomes dense in as h -»• 0. The 
approximate solution, u^ in of (2.1) or (2.3) is defined by 

a(u'^,v'^) = -p„<g,v^>_ for all in S^. (2.5) 

It can generally be proved that u^ u as h 0. For a comprehensive 
treatment of the finite element method, see [10] or [11]. 

Remark 2.3 : It is frequently convenient to replace the equations obtained 
from (2.5) by another system, where all terms that are multiplied by are 
lumped in the diagonal of the matrix. For a discussion of mass lumping, see 
[12]. All numerical results in Section 4 were obtained using a lumped mass 


matrix. 
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Typically, the finite element spaces consist of sufficiently smooth 

piecewise polynomials of some fixed degree defined on a partitioning of D 
into simple subsets with diameter of order 0(h). We have implemented and 
tested a finite element algorithm based on continuous piecewise linear 
functions defined on right triangles. Introducing a basis, bhe 

finite element space (with dimension N = N(h)) in the usual way ([10], 

[11]), the approximate solution u^ in may be expressed as 


u (x,y) = I q^<l»^(x,y). 


j=l 




( 2 . 6 ) 


where = u^(Pj) 
♦j, j=l.‘**,N, we 
column vector g = 


. Substituting (2.6) into (2.5) with given by 

obtain the following system of equations for the unknown 


(q 


r 


■V = 


Aa = s. 


(2.7) 


where g = (g^ 
A is given by 


,gjj)"^ with g^ = and the matrix 


A = (a^j), 'a^j = a(^)^,^j), i,j = 1,‘**,N. (2.8) 

\ 

We close this section by describing a method for efficiently implementing 
the global boundary operator T (consisting of M modes). We start with 
equations (2.7) and let Ny denote the number of grid points in the y 

direction. Note that the iterative solution method to be described in Section 
3 merely requires matrix multiplications of the form Ax with 


T 

X = , so that the matrix A need not be stored. It follows 
readily using an appropriate ordering for the vertices and (2.2), (2.4), and 
(2.8) that for each matrix multiplication Ax we must evaluate 


N 


= I <T(<|> X , 
j=i J 1 J 


1>***»N„ 


(2.9) 


Equation (2.9)' describes a full x Ny matrix with elements 

H,. = <T(iJ) .) ,(}).>„ . If H were an arbitrary full matrix, the work involved 
J J ^ aa 

in computing Ax would be significantly larger than with the local boundary 
conditions, thus degrading the efficiency of the method. However, it is 
apparent from (2.2) that the continuous boundary operator T(u) is of rank 
M and this is also true for the finite element approximation. Therefore 
u = (uj,»»*,Ujj ) can be computed in O(MNy) operations. 

y 

In order to see this, assume for simplicity that M = 1. If we introduce 

T 

the column vector e = (ep»»«,ejj ) defined by e^^ = > then it can 

y CO 


be readily seen that 


Hz = ii, ee z , 


( 2 . 10 ) 


T 

where z = (zj,»»«,Zjj ) is a vector ranging over the boundary points. If 

~ y 

there are M modes in the boundary operator, then Hz is a sum of M terms 


of the form (2.10) 
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3. Solution Method 

As seen In Section 2 , a discretization of (2.1) leads to a system of 
equations 

Ax = b, (3.1) 

where A is typically quite large. Due to the properties of A already 
described, standard iterative methods are not applicable to the solution of 
(3.1), [13]. We have developed an iterative method for the solution of (3.1), 
[1], based on the preconditioned normal equations: 

A'* A'x' = A'*b', (3.2) 

where A' = x' = = Q 2 ^fe» A'* denotes the adjoint of 

A'. The preconditioning matrix M ^ Q 2 ^ was chosen in [1] and [8] to 

be an easily computed partial inverse of Aq, the positive definite matrix 
obtained by setting K = 0 in A. The conjugate gradient method is then 
applied to the system (3*2) and is guaranteed to converge since the normal 
equations are positive definite symmetric# We refer to [1] for a detailed 
description of the conjugate gradient algorithm. The method requires only a 
small number of vector multiplications and additions, and it is only necessary 
to evaluate M~^, A, and A* acting on a vector. Hence no matrices need be 
inverted or stored and the method is efficient provided sufficiently 

reduces the number of iterations. In [1] and [8] we showed this to be the 
case when was chosen to be one or more sweeps of point SSOR or ADI 

applied to the discrete Laplaclan, Aq. 
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Convergence of an iterative method is usually defined by letting the mesh 

size h 0 (hence the number of equations N ®) and keeping all other 

parameters fixed ([13] and [14]), However, in order to maintain accuracy as 

the wave number K increases it is necessary to reduce h. Hence K and h 

must be contrained by means of an accuracy condition, [4], In [4] we showed, 

2 

e.g., that K h must be kept constant for continuous piecewise linear 

finite elements, where in general a > 0 but in many cases a = 0* This 

implies that the number of points/wavelength must be increased as K 

increases* This is also the case as the domain size is increased* 

In view of the above, a more relevant definition of convergence is 
obtained by prescribing a fixed accuracy (e*g*, K^h^ fixed) and letting K 
increase* This has an important bearing on the choice of preconditioner* In 
particular, preconditioners based on fast Laplace solvers are superior with 
respect to the standard definition of convergence but will have a very 

unfavorable growth rate as K Increases* To illustrate this, we consider the 
following model problem on the unit square with Dirichlet boundary conditions 

-Au - K^(l + -i|)u = 0. (3.3) 

The dissipative term 16K models the addition of a term, 6 , to the wave 

equation and is Included to model the radiation condition. 

We will study the convergence rate of the algorithm by analyzing the 

condition number k of the matrix A'*A". It is well known ([14]) that in 

general the convergence rate of the conjugate gradient method applied to (3.2) 

is inversely proportional to A = k ^2 , ^ ^he matrix corresponding 


to the standard five-point discretization of (3.3) (multiplied by h^) 
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Suppose that the preconditioning matrix is positive definite symmetric 

and commutes with A. If = Q2 = M ^^2 , then A is given by 
where is largest eigenvalue of A' (in magnitude) and i® 

smallest eigenvalue of A' (in magnitude). We easily see that the eigenvalues 


X 

max 

^min 


of A are given by 



where 0 < a < 2 and C is a constant. 

a 

'X 0 

We now consider K and h related by K h = e with K h -► 0. The 

eigenvalues are given by 

The cases a = 0 and a = 2 give eigenvalues that are 0(1) and 0(K"’^). 

If we precondition by a complete inverse of the discrete Laplacian, then the 
eigenvalues become 

+~). (3.5) 

The cases o = 0 and a = 2 give eigenvalues that are 0(1) and O(K^). By 
examining the case a = 2/3 we see that the real part of (3.4) and (3.5) can 

vanish. Therefore in (3.4) we can have eigenvalues which are 0(K“^) and A 

o 

can be as large as 0(K ). Similarly, in (3.5) we can have eigenvalues which 

—1 3 

are 0(K ) and A can be as large as 0(K ), thus giving a very slow 


convergence rate for large K. 
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In the remainder of this section we construct a preconditioner based on a 

multigrid method, see (16] , (applied to the discrete Laplacian, Aq) that is 

much more efficient than the preconditioners considered in [1] and [8]. We 

begin by briefly describing a multigrid cycle. Consider a sequence of grids 
0 M 

G , with the mesh size h^ on grid given by 

\ “ 2 ^hQ,i=0,»»* ,M, where h^ = h and hg is independent of h. (Hence 
the number of levels increases as h 0). We obtain a sequence of equations 

= b^, i=l,»»»,M, 


In the same way as (3.1) by discretizing (2.1) as in Section 2 with mesh 
size h^ replacing h. We choose some relaxation scheme for each grid level, 
as well as some interpolation operator from grid G^ to grid G^. 

To begin the cycle, we make r relaxation sweeps on the finest grid 
level, G", and then transfer the residual 


,M-1 


= I 


M-1 

M 


E» . - A« /) 


to grid G^“l . 


On this grid, we obtain the equation 


.M-1 M-1 

A V 


.M-1 


This process is repeated until we get to the coarsest grid G®. On this grid 
we make A relaxation sweeps. To return to the finer grid, we again make % 
relaxation sweeps on G®. We then calculate u^ by 
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1 1 . 0 
U « U , , + I- V . 

new old 0 ~ 


We continue this process, using s relaxation sweeps on each successive grid 
up to the finest grid. The number of relaxation sweeps r,s and I will be 
discussed in the next section. The entire process defines one complete 
multigrid cycle. (See [15] for a more detailed description of the multigrid 
method and [16] for more details on the implementation of multigrid as a 

preconditioner) . 

Since we are using the miltigrid cycle as a preconditioning, we want the 
cycle to be positive definite and symmetric. Assuming a zero initial guess, 
this requires that all operations in the direction of finer grids are the 

adjoint finer grids are the adjoint of the operations in the direction of 
coarsers grids. An efficient relaxation scheme for the Laplace equation is 
Gauss-Seidel with red black (RB) ordering, [9]. This also makes the algorithm 
vectorlzable. In view of symmetry considerations, we see that if RB is used 

when going to coarser grids, then BR must be used when returning to finer 

grids. Furthermore, we maintain symmetry by using linear interpolation for 

and a full weighting for the fine to coarse residual transfer Ij”^» 

[15]. 

It can be shown that a multigrid cycle reduces the error by a fixed 
amount, independent of h, and only requires 0(N) operations. Hence a fixed 
number of multigrid cycles applied to the Laplacian acts as a fast Laplace 
solver. In view of our earlier remark concerning fast Laplace solvers, we 
expect that a multigrid cycle will be effective as a preconditioner for 
small K but the number of Iterations will grow rapidly as K Increases 

subject to the accuracy constraint. These conclusions are confirmed 
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numerically in [16] in connection with a symmetric Helmholtz operator (i.e., a 
Neumann boundary condition replaces the radiation boundary condition). 

In order to obtain a slower growth in the number of iterations as K 
increases, we introduce the following idea. Consider the model problem (3.3) 
and suppose that a preconditioner is used which acts as an inverse of the 
discrete Laplacian (multiplied by h^) on those eigenvectors corresponding to 
a limited part of its spectrum, say those eigenvalues of the form O(h^) with 
0 ^ ^ Oq* Furthermore, assume that the remaining eigenvectors are 

essentially unchanged by the preconditioner, then the eigenvalues of k' are 
given by 

for 0 ^ ^ «Q > 

(3.6) 

for cIq ^ ct ^ 2 
is only 0(K). 

If we now neglect the effects of the residual transfers and coarse to fine 

interpolations in the multigrid preconditioner and assume that the relaxation 

eliminates the highest frequencies of the error on each grid, then we can 

treat multi-grid as a preconditioner of the previous form with eigenvalues of 

the preconditioned matrix given by (3.6). For simplicity assume that 

h * 2“^ on the finest grid. Hence we choose the coarsest grid, G • , to be 

2 ^ 

3 

guch that the highest frequencies on this grid are 0(h ). This gives 



c - CK 2 (1 . 

a K. 




It can be seen that if ~ T * then A 


(3.7) 
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The mimerical examples will demonstrate that for N = 128 three levels are 
optimal. This is vary close to the value predicted by (3.7). Either 
Increasing or decreasing the number of levels results in a striking increase 
in the number of iterations. Furthermore, for N = 200 a simple modification 
of the above argument predicts that a coarsest grid of either 25 or 50 points 
should be optimal. In fact, we find that a coarsest grid of 50 points is 
optimal. Thus, although the above argument is only heuristic, it does 
correctly predict nearly optimal coarsest grids. 


4. Numerical Results 

In this section we present typical numerical results obtained using the 
multigrid preconditioner described in the preceding section. All results were 
obtained for the problem described in Section 2. The results were obtained on 
a square of length u using piecewise linear elements on right triangles. 
There are N grid intervals in each direction, so that the number of 
equations is (N+1)^. In all cases, the mass matrix is lumped. 

The numerical examples are designed to illustrate the convergence 
properties of the algorithm when K increases and K and h are constrained 
by the accuracy requirement, K^h^ fixed. As was indicated in Section 3, for 
large K the preconditioner will be most effective when only a small number 
of levels are used in the multigrid algorithm. The results demonstrate the 
sensitivity of the preconditioner to the number of levels. In addition, the 
numerical results illustrate the effect of boundary conditions on convergence 
properties of the algorithm and the robustness of the method as interfaces 
with strong constrast (causing large backscattering) are introduced. In 
examples 1-4, we assume no interface present and n(x,y) = 1 in (2.1). 
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We define convergence of the iterative method to mean that the normalized 
mean square norm of M”^r is less than 10”^, where r is the residuals 
This quantity is naturally produced by the Implementation of the 
preconditioned conjugate gradient algorithm. We have verified that monitoring 
the norm of r instead of M”^r causes only slight changes in the number of 
iterations. Note that our stopping criteria is more stringent than might be 
required in practice, where the level of the truncation error is used to stop 
the iteration process. 

In example 1 we consider (2.1) with the Neumann data (2. Id) chosen so that 
the exact solution is 


u(x,y) 



(4.1) 


The radiation boundary condition is the local condition, (2.2'), which is 
exact for this mode. Based on the argument of Section 3, we use only three 
grid levels. In addition we use two relaxation sweeps on the coarsest grid 
and one sweep on the other grids. In Table 1, the number of iterations 
required for convergence is shown for different values of K and N. The 
last three entries in Table 1 show the number of iterations required for three 

O O 

frequencies for which K'^h^ = *5 (corresponding to a normalized mean square 
error of about 7%), The first entry corresponds to a cas'e solved using SSOR 
as a preconditioner, which required 284 iterations for convergence, [1]. The 
same problem required 290 iterations to converge using ADI as a preconditioner 
(see [8]). It is apparent that on this simple problem, the multigrid 
preconditioner is more effective than SSOR and ADI. Furthermore, the growth 
in the number of iterations is slow as K increases with K^h^ fixed. 



-18- 


In example 2, we consider the Neumann data as a source modeled by the 

derivative of a Gaussian centered at y “ y • The radiation boundary 

condition is now given by (2.2) with the sum extending over all propagating 
modes and the first four evanescent modes. Results are given for different 
values of K with K^h^ fixed at 1.01 and .425. By examining model 
solutions for these parameters, we believe that the first case correspnds to 
roughly 10% accuracy and the second case corresponds to approximately 5% 
accuracy. The results in Table 2 correspond to three levels in the nultlgrld 
algorithm, with two relaxation sweeps on the coarsest grid and one sweep on 
the other grids. It is seen that the number of iterations grows close to 

O O 

linearly with K for K h fixed. We also observe that for a fixed grid 

size, the number of iterations appears to grow at a rate of 0(K ). This 

shows that the evaluation of the effectiveness of a precohdltloner for large 
K depends crucially on the relationship between K and h. In Table 3, 
results are given using four levels in the multigrid preconditioner and four 
relaxation sweeps on the coarsest grid. It is apparent that the convergence 
is considerable worse in this case. For the last two entries in Table 3, 
seven levels were used and the degradation of the convergence is quite 
striking. 

Due to space limitations, detailed comparisons of the multigrid 

preconditioner with SSOR and other preconditioners will be presented 

elsewhere. We simply state that SSOR for the case of K = 20 and N = 281 

did not converge in 2200 iterations. Furthermore, for K^h^ fixed the number 

of iterations required for convergence is Increasing at a rate greater than 
3 
2 

0(K ). An operation count (counting additions and multiplications equally) 
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Indicates that each Iteration with multlgrld as a preconditioner requires 
about twice the work of an Iteration using SSOR as a preconditioner. The 
0(K) growth In the number of Iterations as K Increases makes the multlgrld 
preconditioner significantly more efficient for large models. 

In the remaining examples we use the same multlgrld preconditioner as In 
example 1. In examples 3 and 4 we consider the effect of the radiation 
boundary condition on the number of Iterations required for convergence. In 
example 3 the Neumann data (2. Id) Is chosen so that the solution Is (4.1), 
while In example 4 the Neumann data Is the derivative of a Gaussian.' In both 
examples K = 16 and N = 201. The radiation boundary condition Is either 
the local boundary condition (2.2') or the global boundary condition (2.2) 
accounting for one, five, ten and twenty modes. The results are presented In 
Tables 4 and 5. In example 4 there are 16 propagating modes, so the only 
boundary condition that accounts for all of the modes Is the last one. In 
example 3, however, there Is only one mode In the solution and all of these 
boundary conditions are non-reflecting on that mode. 

The results In Table 4 Indicate that the radiation boundary condition has 
a very small effect on the number of Iterations required for convergence when 
one mode Is present In the solution. In the case of the derivative of a 
Gaussian, however, the global boundary conditions that allow for reflections 
require considerably more Iterations than the boundary condition that accounts 
for all propagating modes. It can be shown that If the global boundary 
condition (2.2) does not account for all of the propagating modes, the 
boundary value problem (2.1) can be singular or have eigenvalues very close to 
zero. This can degrade the conditioning of the matrix. A comparison of 


Tables 4 and 5 also Indicates that the number of Iterations Increases as the 
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number of propagating modes Increases even if the radiation boundary condition 
is accurate for all of these modes. 

In examples 5 and 6 we consider the effect of a rectangular Interface with 
a piecewise constant index of refraction. In both cases the Neumann data is 
given by the derivative of Gaussian. In example 5, we consider an Interface 
with ^ X Y and 0^7^^ where the index of refraction, n, may be a 
constant other than 1. In example 6 the size of the region is extended, so 
that 7 <x< ^ and 0 ^ y ^ examples n varies from 1 to .25, 
and in one case is 1.25. These constrasts would cause considerable 
backscattering, so that the parabolic equation method is expected to be 
inaccurate for these problems. The number of iterations required for 
convergence is shown in Tables 6 and 7 for these examples. The results 
indicate that the preconditioner is robust and can handle strong constrasts 
extending over relatively large regions. 


5. Conclusions 

We have described a general method to solve Helmholtz type equations for 
an intermediate range of frequencies. The iterative method is based on 
obtaining an effective preconditioner which enables the solution to be 
obtained in a relatively small number of iterations. The relevant limit to 
analyze the convergence properties of a preconditioner is for K increasing 
and a fixed prescribed accuracy level. In this regime, the number of 
iterations increase at a rate greater than O(K^) when using a complete 
multigrid cycle or other methods based on fast solvers as a preconditioner. 
This is vary unsatisfactory for large frequencies. For SSOR, the number of 
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_3 

Iterations increase at a rate greater than O(K^), which is also unfavorable 
for large K, 

The use of multigrid as a preconditioner with a restricted number of 
levels gives a rate of increase of 0(K), thus resulting in a significantly 
more effective algorithm. this is demonstrated by both a heuristic argument 
and by numerical results* Using this method we have been able to solve two- 
dimensional problems with up to 10 wavelengths in each directions and with 
more than 78,000 unknowns in a reasonable number of iterations. For example, 
with a sound speed of 5000 ft. /sec. and a frequency of lOHz, this corresponds 
to a square of length 5000 ft. On an IBM 3033 computer with double precison 
arithmetic, the large model (i.e., k = 20 and 281 x 281 unknowns) required 
four hours of computer time while the smaller model (k = 16 and 201 x 201 
unknowns) required 98 minutes of computer time. The computational effort is 
significantly reduced using the truncation error instead of 10”^ as a 
stopping criterion for the iterative method. Furthermore, both the storage 
requirements and computation time are greatly reduced if we approximate the 
magnitude of the pressure instead of the pressure. This is sufficient for 
many applications. Problems with strong velocity contrast do not appear to 
significantly degrade the performance of the numerical algorithm. Thus this 
method may be suitable for efficient computation of the full acoustic model in 
cases where one-way propagation models and other approximate models would be 
Inaccurate. 
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TABLE 1: Results for Example 1. 


K 

N 


Iterations 

4.16 

61 

.197 

24 

5.92 

65 

.5 

37 

7.76 

97 

.5 

43 

9.4 

129 

.5 

46 



TABLE 2: 

Results for Example 2. 


K 

N 


Iterations 

11.88 

129 

1.01 

640 

16 

201 

1.01 

863 

20 

281 

1.01 

1059 

8.9 

129 

.425 

292 

12 

201 

.425 

472 

15 

281 

.425 

674 
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TABLE 3; Results for Example 2 using 4 levels and 4 relaxation 
sweeps on the coarsest grid. 


K 

N 


Iterations 

11.88 

129 

1.01 

2511 

16 

201 

1.01 

>3900 

20 

281 

1.01 

>2000 

8.9 

129 

.425 

956 

12 

201 

.425 

1492 

15 

281 

.425 

1918 

11.88* 

129 

1.01 

3969 

8.9* 

129 

.425 

1516 

*7 Levels 

and 2 Relaxation Sweeps 

on the Coarsest Grid. 



TABLE 4: Results 

for Example 3 (One Mode 

in Solution). 

Boundary Condition 

Modes in BC 

Iterations 

Local 

— 

81 

Global 

1 

77 

Global 

5 

80 


Global 10 81 

Global 20 81 
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TABLB 5: Results for Example 4 (Point Source). 


Boundary Condition 

Modes In BC 

Iterations 

Local 

— 

805 

Global 

1 

1401 

Global 

5 

1411 

Global 

10 

1226 

Global 

20 

863 




TABLE 6: 

Results for Example 5* 


K 


N 

Index of 
Refraction 

Iterations 

16 


201 

1 

863 

16 


201 

0.5 

856 

16 


201 

0.33 

798 

16 


201 

0.25 

841 
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TABLE 7: Results for Example 6 


K 

N 

Index of 
Refraction 

Iterations 

16 

201 

1 

863 

16 

201 

0.5 

943 

16 

201 

0.33 

1108 

16 

201 

0.25 

788 

16 

201 

1.25 

1082 

20 

281 

1 

1059 

20 

281 

0.33 

1269 
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